require(runjags)
require(coda)
require(stringi)

source("A10_functions.r")

data <- read.csv("C1_main_data.csv")
data$date <- as.Date(data$date)

#### estimation ####
# preprocessing
n <- nrow(data) - 1
m <- 4
y <- cbind(data$L_K_3s, data$L_D_3s, data$L_C_3s, data$L_L_3s)[-1, ]
k <- rowSums(y)
y[k == 0, ] <- NA
k[k == 0] <- 1

x <- cbind(diff(data$cabinet / 100), 
           diff(data$LDP / 100), 
           data$DPJ.gov[-1], 
           diff(data$cabinet / 100) * data$DPJ.gov[-1], 
           data$HOR[-1] / 47, 
           (data$HOR[-1] / 47) ^ 2, 
           data$HOC[-1] / 35, 
           (data$HOC[-1] / 35) ^ 2, 
           diff(data$pmid))
l <- ncol(x)

LDP.mean.minus.SD <- mean(x[, 2]) - sd(x[, 2])
LDP.mean.plus.SD <- mean(x[, 2]) + sd(x[, 2])

# initial values
initial.values <- list()
for (i in 1:3) {
  initial.values[[i]] <- inits.fun(10 * i, 100 * i, i)
}

# estimation
result <- run.jags(model = "B1_dynamic_multinomial_model.R", 
                   monitor = c("beta", "Omega", "theta", "theta.star"), 
                   data = list(y = y, k = k, x = x, n = n, m = m, l = l, 
                               mu0 = rep(0, m - 1), 
                               Omega.star0 = diag(m - 1) * 0.01), 
                   inits = initial.values, modules = "glm", 
                   n.chains = 3, burnin = 5000, sample = 5000, 
                   adapt = 20000, summarise = FALSE, 
                   thin = 50, method = "parallel")
# save(result, file = "D6_LDP_result.Rdata")
# load("D6_LDP_result.Rdata")

# convergence check
sum(gelman.diag(result$mcmc, multivariate = FALSE)$psrf[, 1] > 1.1)

# combine MCMC draws into a single matrix
sims.matrix <- rbind(result$mcmc[[1]], result$mcmc[[2]], result$mcmc[[3]])
colnames(sims.matrix) <- colnames(result$mcmc[[1]])

#### results ####
var.names <- c("Change in cabinet approval", "Change in LDP support", 
               "DPJ government", "Cabinet x DPJ gov.", 
               "Time since HoR election", 
               "Time since HoR election squared", "Time since HoC election", 
               "Time since HoC election squared", "Inauguration of the PM")

# Table A.8
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], 
                        var.names), 3)
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], 
                        var.names), 3)
round(posterior.summary(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], 
                        var.names), 3)

#### simulation for the LDP support rate ####
LDP.seq <- seq(round(min(x[, 2]), 3), round(max(x[, 2]), 3), 0.001)

# Figure A.6
sim.result.LDP <- array(NA, c(length(LDP.seq), 4, 3))
for (i in 1:length(LDP.seq)) {
  sim.x <- x
  sim.x[, 2] <- LDP.seq[i]
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1 <- exp.theta.star.1 / sum.exp.theta.star
  prob.2 <- exp.theta.star.2 / sum.exp.theta.star
  prob.3 <- exp.theta.star.3 / sum.exp.theta.star
  prob.4 <- 1 / sum.exp.theta.star
  sim.result.LDP[i, 1, ] <- posterior.summary(rowMeans(prob.1))
  sim.result.LDP[i, 2, ] <- posterior.summary(rowMeans(prob.2))
  sim.result.LDP[i, 3, ] <- posterior.summary(rowMeans(prob.3))
  sim.result.LDP[i, 4, ] <- posterior.summary(rowMeans(prob.4))
}

png("Figure_A6.png", 
    width = 6, height = 1.6, units = "in", pointsize = 7, res = 1500)
layout(matrix(1:4, 1, 4))
par(mar = c(4, 4, 2, 1), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0.1, 0.3), 
     main = "Komeito", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 1, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0.1, 0.3), 
     main = "DPJ", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 2, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0, 0.2), 
     main = "JCP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 3, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.053, 0.089), ylim = c(0.45, 0.65), 
     main = "LDP", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = LDP.mean.minus.SD, lty = 3, col = "gray80")
abline(v = LDP.mean.plus.SD, lty = 3, col = "gray80")
line.fun(LDP.seq, sim.result.LDP[, 4, ])
rug(x[, 2])
box()
mtext("Change in LDP support", 1, 2.5, cex = 0.7)
mtext("Classification probability", 2, 2.5, cex = 0.7)
axis(1, at = seq(-0.04, 0.08, 0.02), 
     labels = c("-0.04", "", "0.00", "", "0.04", "", "0.08"), lwd = 0.5)
axis(2, lwd = 0.5)
dev.off()

# difference in predicted probabilities between cases with changes in the LDP public support rate one SD above vs. below its mean
sim.diff.LDP <- function(low.value, high.value, prob = 0.95) {
  sim.x <- x
  sim.x[, 2] <- low.value
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.low <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.low <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.low <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.low <- 1 / sum.exp.theta.star
  sim.x[, 2] <- high.value
  exp.theta.star.1 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),1\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[1,")], sim.x))
  exp.theta.star.2 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),2\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[2,")], sim.x))
  exp.theta.star.3 <- exp(tcrossprod(rowMeans(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "star\\[([0-9]|[1-9][0-9]|[1-2][0-9][0-9]),3\\]")]), rep(1, nrow(sim.x))) + 
                            tcrossprod(sims.matrix[, stri_detect_regex(colnames(sims.matrix), "beta\\[3,")], sim.x))
  sum.exp.theta.star <- exp.theta.star.1 + exp.theta.star.2 + exp.theta.star.3 + 1
  prob.1.high <- exp.theta.star.1 / sum.exp.theta.star
  prob.2.high <- exp.theta.star.2 / sum.exp.theta.star
  prob.3.high <- exp.theta.star.3 / sum.exp.theta.star
  prob.4.high <- 1 / sum.exp.theta.star
  prob.1.diff <- prob.1.high - prob.1.low
  prob.2.diff <- prob.2.high - prob.2.low
  prob.3.diff <- prob.3.high - prob.3.low
  prob.4.diff <- prob.4.high - prob.4.low
  sim.diff <- rbind(posterior.summary(rowMeans(prob.1.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.2.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.3.diff), prob = prob), 
                    posterior.summary(rowMeans(prob.4.diff), prob = prob))
  rownames(sim.diff) <- c("Komeito", "DPJ", "JCP", "LDP")
  colnames(sim.diff) <- c("Est", "Upper", "Lower")
  sim.diff
}

round(sim.diff.LDP(LDP.mean.minus.SD, LDP.mean.plus.SD), 3)